Heated granular fluids: the random restitution coefficient approach 
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We introduce the model of inelastic hard spheres with random restitution coefficient a, in order 
to account for the fact that, in a vertically shaken granular system interacting elastically with the 
^3 , vibrating boundary, the energy injected vertically is transferred to the horizontal degrees of freedom 

through collisions only, which leads to heating through collisions, i.e. to inelastic horizontal collisions 
with an effective restitution coefficient that can be larger than 1. This allows the system to reach 
, a non-equilibrium steady state, where we focus in particular on the single particle velocity distri- 

bution f(v) in the horizontal plane, and on its deviation from a Maxwellian. Molecular Dynamics 
simulations and Direct Simulation Monte Carlo (DSMC) show that, depending on the distribution 
of a, different shapes of f(v) can be obtained, with very different high energy tails. Moreover, the 
fourth cumulant of the velocity distribution quantifying the deviations from Gaussian statistics is 
obtained analytically from the Boltzmann equation and successfully tested against the simulations. 
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Despite a wealth of recent experimental investigations, it seems that there is no consensus yet as to the characteristics 
of the collective behaviour of vibrated granular mono- or multi-layers . Among the statistical properties of interest, 
the velocity fluctuations and more precisely the deviations from Maxwell-Boltzmann distribution due to inelastic 
collisions have been intensively studied for rapid granular flows, and several groups reported an overpopulated high 
q ■ energy tail that can be fitted by stretched exponentials with various exponents. 

Theoretical investigations commonly use the inelastic hard spheres (IHS) model which has proven very useful 
■ despite the simplicity of its definition: smooth hard spheres undergo binary inelastic momentum conserving collisions, 
thereby losing during the collision a constant fraction of their relative normal velocity (and therefore loosing energy) . 
Many studies have concentrated on the homogeneous cooling state of the IHS [^0, obtained by letting the system 
I/"") | evolve without any energy injection. On the other hand, experiments on strongly vibrated granular media consider 
systems that are heated by an external forcing (e.g. a vertically oscillating plate) to compensate for the energy loss 
due to collisions, and are therefore in a non equilibrium stationary state (NESS). Various ways of modeling the heating 
mechanism have been put forward, mostly consisting in the study of an effective system with a constant coefficient 
| of normal restitution a, in which the energy lost through collisions is balanced by energy injected by random "kicks" 
|^5|-|l9|. In this context, analytical and numerical results have been obtained for the high energy tail of the velo city 
distributions and for its fourth cumulant, a natural measure of the deviation from a Gaussian distribution fl^ , [l7| , po| . 
Moreover, more realistic Molecular Dynamics simulations of vibrated soft sphere mono-layers have been able to 
reproduce the phenomenology obtained in some experiments ||. 

In this paper, we propose an alternative modeling of a three-dimensional system shaken vertically along the z axis, 
for which the velocities are studied both numerically and analytically in the horizontal (xy) plane. If the collisions 
with the boundaries are elastic, the vibrating wall feeds energy in the system in the z direction only, but not in the 
xy plane where the correponding two dimensional energy is either lost or gained at each collision between two grains. 
Upon restricting to the xy velocities, we obtain a system that is subject to an effective stationary dynamics with 
sequential energy injection or dissipation consecutive to particle-particle collisions only. We shall thus disregard the 
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■ vibrating boundary and concentrate on the effective planar dynamics where energy gains and losses statistically cancel 
in the NESS, so that an effective restitution coefficient can either be larger or smaller than 1; this will be accounted for 
by a random a, drawn from a probability distribution p(ct) (such that the second moment a 2 = 1 in order to conserve 
energy on average). Our approach, in which momentum transfer occurs only through collisions, consequently differs 
from the situations investigated in Jl5|-p^| where energy is injected globally into the system at regular time intervals 
through a stochastic external force. Clearly, the functional form of p(a) reflects the energy injection mechanism, and 
it is a difficult task to establish this connection. The physical situations that could be described by the present model 
are mono- or multi-layers of grains provided the interactions with the wall are elastic, or a multilayer system with 
arbitrary interactions with the wall, in a "bulk" region far from the boundaries where energy and density can be 
considered as constant. In a very different context, a similar stochastic coefficient of restitution has been introduced 
to study the dynamics of a one-dimensional granular gas, thereby accounting for internal degrees of freedom |22|. 
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The paper is organized as follows: after setting the general framework in section II, we consider the two-dimensional 
projection of three-dimensional simulations of inelastic hard spheres with constant restitution coefficient a, energy 
being injected in the third direction (section III). We analyze in particular the distribution of energy transfer through 
collisions. In section IV, the two dimensional IHS model with random restitution coefficient is then studied by 
analytical and numerical means. Section V is devoted to a short investigation of the one-dimensional case, and some 
conclusions are finally presented in section VI. 



II. GENERAL FRAMEWORK 



A. IHS: definitions and notations 



In the inelastic hard spheres (IHS) model, grains are modeled as smooth hard spheres of mass m undergoing binary, 
inelastic and momentum-conserving collisions: a collision between two spheres labeled by 1 and 2, with velocities vi 
and V2, dissipates a fraction (1 — a) of the component of the relative velocity V12 = vj — V2 along the center-to-center 
direction <x. Noting with stars the post-collision velocities, this translates into v* 2 <r = —a Vx2 <r, while the tangential 
relative velocity (perpendicular to a) is conserved, i.e.: 

vj = vi - - (1 + a) (vi2 • <t)ct 

v 2 = v 2 + - (1 + a) (vx2 • <t)ct . (1) 

These equations are also sometimes written in terms of restituting collisions, i.e. for collisions which yield (vx, v 2 ) as 
postcollisional velocities, for precollisional velocities (v**, v"): 

v** = V! - i ^1 + (vx2 • a)a 

vj* = v 2 + \ (l + ^ J (v 12 • <x)«r . (2) 



B. Evolution equation 



The Enskog-Boltzmann equation describes the evolution of the one-particle distribution function /(r, v, t), upon 
the molecular chaos hypothesis p3fl . In the homogeneous case, for hard spheres of diameter a, in d dimensions, this 
equation reads: 



df(vi,t) 



X(J d 1 J ' dv 2 J d&(vi 2 ■&) 



/(vr,t)/(vr,o-/(vi,t)/(v 2 ,t) 



at J Z J v 1Z ' l« 2 ' 

= Xl(f,f) (3) 

The prime on the integration symbol is a shortcut for J d&Q(vi2 ■ <r) where 6 is the Heavyside function, while x 
accounts for excluded volume effects (for elastic hard spheres, \ coincides with the density dependent pair correlation 
function at contact). The forcing mechanism necessary to sustain a NESS can be of different types, and has not been 
introduced in Eq. Q3[) . This issue will be addressed in section IV. 

C. Sonine expansion for the velocity distribution 

Because of analytical, numerical and experimental evidences, it is customary to look for scaling solutions of eq. (|^), 
in the form |26| 

/( v ^) = ^7Tt/(^77t) . ( 4 ) 
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where n is the density and the thermal velocity vo is by definition related to the temperature T(t) through -?t>o(i) = 
T(t), where in turn the temperature is defined by the average kinetic energy of the particles: 

Replacing the scaling function in eq (|^), together with the law of evolution of the temperature (see for details) 
yields the following equation: 

111 ^d + c 1 ^-)f(c 1 )=I(f~,f) , (6) 



d \ dci 
where c< = Vi/v (t), 

hi f) = J dc 2 J' da(c 12 ■ a) I -L/ (cr) ; (cr) - f( Cl )f(c 2 )\ , (7) 

and n p = -f dc 1 c^/(/, /). 

In the elastic case (a = 1), / is the Gaussian $(c) = ■K~ d l' 1 cxp(— c 2 ) (since /(<&,$) = 0). The deviations from <!> 
are studied by an expansion in terms of Sonine polynomials [^4| 

/>)=*(c)(l + f>A(c 2 )^ . (8) 

Indeed, the coefficients a p can be obtained from the moments of /, because the S p satisfy the orthogonality relations: 

dc<S>{c)S p {c 2 )S pl {c 2 ) = 5 PP ,N P (9) 

where the N p are normalization constants. These polynomials consequently depend on the dimension d; the first ones 
read: 

So (a?) = 1 
Si(x) =- x +2 

„ , x x 2 d + 2 d(d + 2) 

s ^ = Y~^r x+ ii (10) 



The orthogonality properties allow to write 



a P = j^(S p (c 2 )). (11) 



In particular, ai = (2/d)(—(c 2 ) + d/2) vanishes because of the definition of temperature, and a 2 is related to the 
fourth cumulant of /: 

4(c 4 ) 

"'W^- 1 - (12) 

In section IV, we shall briefly recall the different steps involved in the computation of a 2 that quantifies the deviations 
from gaussianity. 

D. DSMC Method 

The Direct Simulation Monte Carlo (DSMC) method allows to solve numerically the Boltzmann equation (||) p7[ . 
It has been successfully used e.g. for the study of the homogeneous cooling state of the IHS, in |ll]] to assess the 
validity of the Sonine expansion, and in j2f| to analyze the high energy tail of /. In p3), the case of heated granular 
fluids has also been considered and compared to the theoretical predictions of |l7]] . 
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The positions of the particles do not appear in equation (|3j). Therefore all reference to space is here useless; this 
amounts, in the usual DSMC language, to taking only one cell where all particles stand, and of course eliminates the 
possibility to study spatial inhomogeneities. Moreover, upon an appropriate rescaling of time, we shall take % = 1 
and cr = l. Since no external force is present, the velocities of the particles do not change between collisions, and 
no "free streaming stage" is here necessary, and the simulation has the following scheme: we start from N particles 
with random velocities taken from an arbitrary distribution (e.g. Gaussian, or flat); then the evolution proceeds by 
choosing at random pairs of particles a direction <x for their center to center direction, and updating their 

velocities according to the collision rule with a probability proportional to 6(v^ • <x)vjj • a. Once a stationary state 
is reached, running averages can be taken on the quantities of interest. 



E. Molecular Dynamics simulations 



The Molecular Dynamics (MD) simulations integrate the exact equation of motion of the model, with no reference 
to the Boltzmann equation: we consider N spheres of diameter er, in a box of linear size L in dimension d, with 
periodic boundary conditions |p8| , p9| . These spheres initially have random velocities, and we use an event-driven 
algorithm to study their dynamics. Once more, running averages are taken once a stationary state is reached. 



III. 2D PROJECTION OF A 3D SYSTEM 



We consider an IHS model with constant (velocity independent) restitution coefficient a < 1, in dimension d = 3, 
the Cartesian coordinates being labeled x,y and z. The collisions are dissipative but energy is re-injected in the 
vertical direction (z) in the following way: after each collision, the z-component of the velocity of the particles having 
collided is randomly drawn from a Gaussian distribution at fixed temperature, while the x — y components remain 
those resulting from the initial inelastic collision. This procedure is intended to mimic both the energy injection due 
to a vibrating wall, and transfer to horizontal degrees of freedom through collisions only: indeed, a vibrating boundary 
can yield in the bulk of a multi-layer system an equilibrated Gaussian vertical velocity ]2|,[50| . We could have used 
different vertical velocity distributions, in particular asymetric ones, with the same conclusion: the point here is to 
illustrate the horizontal energy tranfer mechanism. 

Let us indeed analyze the system in the horizontal plane (xy). After a transient, the two-dimensional temperature 
T X y — \{{ v lc JrV y)) remains stationary (see the inset of Fig. |l|): overall, i.e. when the three components of the velocities 
are considered, a collision is dissipative, but in the rry-plane, energy can be gained by a transfer from the ^-direction. 
In the horizontal plane, both phenomena compensate each other, as appears in Fig. 0: the histogram of the energy 
transferred in the a;y-plane at each collision (given by AE = v'i x + v'i y + v'£ x + v'£ y — (v\ x + v\ y + v\ x + v\ y ), where 
1 and 2 label the two colliding particles) has of course a negative part but also a positive one, corresponding to an 
energy gain in the xy-plane thought a three-dimensional collision (the positive part would be absent in a simulation 
without heating and constant a < 1). This positive part compensates the energy loss and allows to reach a NESS 
with a constant T xv . 
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FIG. 1. Histogram of the two-dimensional energy change consecutive to a collision for a 3-D DSMC simulation of a system 
of TV = 5 10 5 particles, with constant restitution coefficient a = 0.9 and energy injection in the z-direction. The positive 
part corresponds to effective restitution coefficients larger than 1. The inset shows the transient approach to the NESS of the 
two-dimensional temperature for a = I, 0.9, 0.6 (from top to bottom). The initial value of the temperature coincides with 
that of the equilibrated vertical degrees of freedom. A similar anisotropy between horizontal and vertical temperatures has 
been observed in |H| 

It would be feasible to inject energy in a more refined way, by considering for instance the collisions with the 
vibrating wall. Seeking for analytical results concerning the 2D horizontal velocity distribution, we shall however make 
the simplifying assumption that the effective 2D restitution coefficient decouples from the impact relative velocity. 
The resulting zeroth order modeling introduced below displays qualitatively the same energy transfer behaviour as 
the projected 3D system and is amenable to a kinetic theory description. 



IV. RANDOM RESTITUTION COEFFICIENT 



We shall hereafter consider an IHS model for which, at each collision, the restitution coefficient is drawn from a 
probability distribution p{a). The means over p(a) will be denoted by an overline, and the distribution of a 2 by 
p(a 2 ). Since, in a binary collision with restitution coefficient a, the energy change is: 

A£=f(v^) 2 ^, (13) 

(where is the relative velocity before collision, & the center to center direction, and m the mass of the particles), 
and since the value of a is taken uncorrelated with the velocities of the particles, we shall consider distributions 
with a 2 — 1 in order to ensure a stationary, constant temperature regime (at each collision, energy changes, but is 
conserved on average: AE = 0). Moreover, we will restrict ourselves to positive values of a. Since the average energy 
is constant, the granular temperature is also a constant determined by the initial velocity distribution. This model is 
therefore intended to study the distribution of rescaled velocities c ~ v/«o- 



A. Analytical results 

The methods of (Tt]] for the case of constant normal restitution (with or without external heating) can be easily 
applied to the case of random a to systematically obtain the coefficients of the Sonine expansion. We refer to |l7]] 
for details and sketch only the principal steps. Note that the calculation is made under the hypothesis of a small a 2 : 
the Sonine expansion (||) is therefore truncated after second order and besides, terms of order a 2 are discarded (it is 
possible to go beyond the linear approximation in with again truncation of (j^) after n — 2; it was however shown 
that in the case of constant a, the correction is less than 10% fl3|). 

Once the moments p p = — J dci c\I(f, /) have been denned, multiplying the equation (|^) by c\ and integrating 
over ci yields 

P P = P2- d {c p )- (14) 

Taking p — 4 and approximating / by its second order Sonine expansion, it is now possible to evaluate /i2, (14 and 
(c 4 ), that can be averaged over a: 

(c < ) = fS±2 (1 + <l2) 

JM = \f^n d (7\ + a 2 n) , (15) 

with Tx = (1 - a 2 )(d + 3/2 + a 2 )/4, T 2 = 3(1 - a 2 )(10d + 39 + 10a 2 )/128+ (1 + a)(d - l)/4 (Q d is the volume of 
the d-dimensional unit sphere). Inserting these relations in (14) and neglecting terms of order a 2 , leads to the final 
expression for 02 in dimension d: 
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1fi 1 - 3a 2 + 2a 4 

a2 = 16 = = . (16 

9 + 24d + 32(d - l)a + (8d - ll)a 2 - 30a 4 

It can be checked that the expression obtained in [ p7j is recovered in the case of constant a: when = £(a — a*), 
we get 

= lf 1 - a* - 2a 2 - + , , 

° 2 9 + 24d + 8o!*d-4lQ!* + 30(l-a*)a2 - 1 > 



Moreover, in the case a 2 = 1, expression (16) reduces to 

a 4 - 1 

a 2 = 16 = . (18) 

16d-l + 16(d-l)a-15a 4 

The values of a 2 corresponding to the various distributions of normal restitutions shown in Fig. ^, are given in 
appendix A. 

An important consequence of Eq. (|l6|) is that the fourth cumulant depends only on a, a 2 , and a 4 : two different 
distributions pip) having the same first, second and fourth moments should then yield very similar velocity distribu- 
tions. This will be checked numerically in the next subsection and can be considered as a test for the consistency of 
the linear order approximation in a 2 underlying the analytical computation. 
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FIG. 2. Various distributions p(a) of the restitution parameter, fulfilling a 2 — 1: 

a) trimodal distribution p(a) = §[5(a — y/1 + 7) + S(a — y/l — 7)] + (1 — b)8(a — 1), with the bimodal case included for 6=1; 

b) flat p(a) (a G [0.457427; 1.457427]); c) flat p(a 2 ) (i.e. linear p(a)), for a 2 G [1 - 7, 1 + 7]; d) Gaussian p(a 2 ). 

The question of the high energy tail, which was addressed in (l^Jl^l by neglecting the gain term in the collision 
integral /, turns out to be more problematic. Whenever a is allowed to take values exceeding 1, the gain term can 
no longer be discarded and we could not obtain analytical predictions. We shall therefore resort to a numerical 
resolution of the Boltzmann equation and to molecular dynamics simulations to analyze the high velocity statistics. 
It can however be shown that the velocity distribution is Gaussian in the elastic case only [i.e. for p(a) = 8(a — 1), 
see appendix B] . 



B. Numerics 



We have simulated the IHS model with random restitution coefficient in two dimensions for various distributions 
p(a), both with DSMC and MD methods. The DSMC simulations were performed with 3.10 5 particles, the MD ones 
with up to 50000 particles, and packing fractions from 10% to 40%. The restitution coefficient a were drawn from 
distributions of various types: bimodal or trimodal distributions, flat distributions, flat distributions of a 2 , Gaussian 
distributions of a 2 (with a cutoff in zero to ensure that a 2 > 0). The distributions are shown in Fig. ^: no large 
values of a or pathological distributions will be used. 
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The stationarity of the kinetic energy (following from a 2 — 1) is controlled during the simulation, and allows to 
obtain in a single run the velocity distribution with a precision of typically 7 to 8 orders of magnitude. We first show in 
Fig. H the histogram of the energy transfers during the collisions, corresponding to various p(a), for both DSMC and 
MD simulations; the similarity with the results of the 3 dimensional simulations (Fig. allows to validate the model 
as far as energy transfer is concerned (the precise form of the histogram depends on the probability distributions p(a), 
but the qualitative features are those displayed by the projected 3D system considered in section III). 
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FIG. 3. MD (symbols) and DSMC (lines) simulations of the random a model. MD curves with packing fraction 
10%. The various curves correspond to various p(a): from top to bottom (on the right), flat p(a) (DSMC), bimodal 
p(a 2 ) = \{8(a 2 - 1/2) + S(a 2 - 3/2)) (DSMC and MD), uniform a 2 G [0.5; 1.5] (DSMC and MD), and bimodal 
p(a) = |(5(a - 1.04) + 5 {a - 0.958332) (DSMC). 

Figure || shows the velocity distributions for a flat p(a 2 ) between and 2, for DSMC and MD simulations. The 
agreement between both sets of data is remarkable, and was also checked for other choices of p(a). Moreover, the 
curves obtained in MD simulations with small or large packing fractions (up to 40%) are indistinguishable (not shown). 
The inset shows that the distribution of impact parameters in Molecular Dynamics simulations is flat, which is a hint 
that no violation of molecular chaos is observed and an indication that the factorization of the 2-particle correlation 
function resulting in Eq. ([}]) holds. The super-imposition of MD (where a priori inhomogeneities and/or violations 
of molecular chaos could appear), and DSMC results is in contrast with the phenomenology at constant a (3l]] or 
with randomly driven IHS |32]| , and is probably due to an efficient randomization of the velocities with the collision 
rule of the present model. For a constant dissipative restitution parameter, colliding particles emerge with more 
parallel velocities than in the elastic case a = 1. When they recollide, their velocities are still more parallel. Here, 
this mechanism for the creation of velocity correlations violating molecular chaos seems removed by the possibility of 
having a > 1. This validates the theoretical approach based on the Boltzmann equation. 

The remainder of this article is devoted to the deviations from gaussianity that can be seen in Fig. f§ To this 
aim, we shall use DSMC simulations that allow to obtain precise velocity distributions for a larger range of velocities 
than the MD method. However, as stated above, we always observed an excellent agreement between DSMC and MD 
velocity statistics, up to the resolution of MD. 
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FIG. 4. Velocity distributions for the random a model, for MD (circles) and DSMC (solid line) simulations with respectively 
50000 and 300000 particles, with the same /9(a) [flat p~(a 2 ) between and 2]. The agreement between both sets of data is 
striking. Large deviations from the Maxwellian (dotted line) are observed. The inset is the distribution of impact parameters 
for the MD simulations, showing no violation of molecular chaos. 

The theoretical (from eq (|l6|)) and measured (from the fourth moment of the velocity distribution) values of a 2 are 
given in Table 1. For small a 2 the measured and theoretical values agree perfectly, while for larger 02 the agreement is 
worse, probably due to the neglect of quadratic a\ and higher order terms in the derivation of the theoretical formula. 
In fact the disagreement is at most of 10% when \d2\ > 0.1, comparably to the results for constant a. As predicted, 
the same value of the fourth cumulant (12 is associated with distributions having the same first, second and fourth 
moments [see Eq. (|l8|)]. Moreover, Fig. || shows strong numerical evidence that, in this case, the whole distributions 
of velocities are indistinguishable, i.e. not only ai but the complete statistics (including large velocities) depend only 
on the first moments of p(pi). Of course, for other p(a), this property may not be true anymore. 

In fig. we show the comparison between the normalized velocity distribution for two distributions of restitution 
coefficients, together with the Sonine prediction 1 + a^Si- In the case of small 02, the agreement is striking: the 
Sonine expansion is supposed to be valid for small values of c, while we see in Fig. [| an agreement up to relatively 
large c. For larger 02, the global shape corresponds to the prediction, but the quantitative agreement is lost. This is 
quite expected since, when ai is not very small, higher order terms in the Sonine expansion become relevant. 




FIG. 5. Velocity distribution functions for DSMC simulations of the random a model, for two sets of two distribu- 
tions having the same values of a and a 4 : Trimodal p(a) = §(5(a — \/l + 7) + 8{ol — ^1 — 7)) + (1 — b)S(a — 1) with 
(b ~ 0.546248, 7 « 0.390584) (solid line) versus a flat distribution for a 2 € [0.5; 1.5] (circles), and another trimodal distribution 
with (6 « 0.47779,7 ~ 0.835254) (dotted line) versus another flat distribution for a 2 £ [0 : 2] (squares). The agreement over 8 
orders of magnitude shows that the 02 approximation yields reliable predictions. 



8 




0.8 



FIG. 6. Comparison of the ratios //$, measured in DSMC simulations (symbols), versus the theoretical predictions of the 
Sonine approximation 1 + 0,282 (lines), for two flat distributions of a 2 : a 2 £ [0,2] (circles) yields a relatively large fourth 
cumulant (02 ~ 0.18) and the agreement is satisfying at small velocities only, while the measured and theoretical curves are 
indistinguishable for a 2 £ [0.5; 1.5] (crosses), for which 02 ~ 0.04 is small. 

We now turn to the study of the large velocity tails. A first indication is given by the plot of the derivative of 
In /(c), which is linear for a Gaussian, and constant for an exponential law. An exp(— Ac B ) tail on the other hand 
leads to a c B_1 behaviour. We see in Fig. that the non gaussianity is indeed revealed by this criterion, but that 
the numerical noise hinders any clear conclusion on the value of B. Since our velocity statistics are smooth over 8 
orders of magnitude, this approach is unable to determine the values of B, even if f(v) behaves asymptotically as a 
stretched exponential. 

Figure || on the other hand shows three fits, for three distributions of restitution coefficients. These fits are of the 
form exp(— Ac B ), and we obtain a wide range of possible values for B Q from 0.8 to 2, with fits accurate over 6 
orders in magnitude. In particular, a convenient choice of p(a) is compatible with B = 1.6, which has been found in 
some experiments (close to B = 3/2 obtained in [^7) for randomly driven IHS fluids). The corresponding p(a) 
[trimodal p(a) = \{8{a--s/T+^)+8{a--s/T^)) + {l-b)8{a-l) with (b « 0.546248, 7 w 0.390584), or equivalently 
flat a 2 € [0.5, 1.5]], is not very broad and does not imply the use of particularly large values of a. Other choices of 
p(a) can also lead to an exponential tail (B = 1), similarly to the case of the homogeneous cooling state (constant 
a < 1, with no heating), or even to larger distributions with B < 1. 
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1 As stated above, for the situations investigated, the value of B seems to depend on the first moments of p{a) only. 
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FIG. 7. Derivative d\nf/dc versus c for various p(a): from top to bottom on the right, flat p(a 2 ) with a 2 £ [0; 2], bimodal 
distribution p(a 2 ) = \ [5{a 2 - 1/2) + S(a 2 - 3/2)], and fiat p(a 2 ) in [0.5, 1.5]. 















































CIO 


nnnnmriin in 







FIG. 8. Velocity distributions for the same p(a) as in Fig. g (circles: flat p(a 2 ) in [0.5,1.5]; diamonds: bimodal a 2 — 1/2 
or 3/2; squares: flat p(ct 2 ) in [0; 2]), obtained from DSMC simulations, together with the Maxwellian <E> (solid line) and fits to 
K exp(— Ac B ) (dotted line: exp( — 1.5c 1 ' 6 ); dashed line: 10 exp(— 3.4c); dot-dashed line: 100 exp(— 5c 0,8 )). The fits are accurate 
over 6 orders of magnitude and values of B consistent with the experimental data (B ~ 1.6) can be obtained for a convenient 
choice of p(a). 

V. ONE-DIMENSIONAL CASE 

For completeness, we briefly report the results obtained for the one-dimensional version of the random-a model; 
it should correspond to a projection in one dimension of a two-dimensional system, with energy injection along the 
projection direction. 

The study of the projected model defined in section III (injection of energy by randomly drawing the y-component 
of the velocity after each collision) yields similar results as those obtained in the case of a three-dimensional system 
projected in 2D. 

The random-ct model however displays a certain number of pathologies: for d = 1 and a 2 = 1, we obtain analytically 
a2 = —16/15 from Eq. (|l8|) . This large value, independent of p(a) (i.e. non perturbative), indicates that the Sonine 
expansion is not valid. In fact, numerical investigations (both MD and DSMC) show that the velocity distributions 
f(v) have power law tails (see Fig. [)]). This behaviour differs significantly from a Gaussian and explains the failure 
of the Sonine approximation. 
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FIG. 9. Velocity distributions for the one-dimensional case, for various p(a): fiat p (middle), flat p (a 2 £ [0.5; 1.5]) (top), 
and bimodal p(a) = \{8(a - 1.04) + 5 (a - 0.958332) (bottom). Data for DSMC (symbols) and MD (lines) simulations are 
shown. Some distributions have been shifted for clarity. 



VI. SUMMARY AND CONCLUSIONS 



We have introduced the idea of a random restitution coefficient in the IHS model to describe a vertically vibrated 
layer of granular material. Since energy is injected only along the vertical axis, and transferred through collisions in 
the perpendicular directions, our approach accounts for the fact (as shown in section III) that the projection in 2 
dimensions of a 3-dimensional collision can correspond to a gain in the two-dimensional energy, and therefore to an 
effective restitution coefficient a larger than 1, even if the genuine a necessarily corresponds to a dissipative collision. 
The model is consequently studied in 2 dimensions, with a probability distribution p{a) for the restitution coefficient. 
We have analyzed the velocity distributions, and in particular the deviation from the Maxwellian, using the Soninc 
expansion technique. Following JL7| we obtained analytically the expression of the fourth cumulant which has 
been tested against Molecular Dynamics (MD) and Monte Carlo Direct Simulations (DSMC). It turns out that the 
theoretical predictions for ai are quite accurate, with a slight overestimation for 0,2 that probably corresponds to the 
approximations made during the calculation (nonlinear terms 0{a\) and higher order Sonine polynomials neglected); 
in particular a 2 depends only on the first moments of p(a); numerically, the whole velocity distribution has the same 
property with a very high precision, at least for the p{a) studied here (we do not exclude that this behaviour could 
be violated by other kinds of distributions). Moreover, the comparison between numerical data and the second order 
Sonine expansion shows a remarkable agreement for small values of ai- The high energy tails, studied with DSMC 
simulations, can be fitted by functions of the form exp(— Ac B ), with B < 2 depending on p(pt). It would certainly 
be interesting to have theoretical predictions concerning B\ it seems for example that the high-energy tail is always 
overpopulated with respect to the Maxwellian, while a particular choice of heating can also yield an under-population, 
as shown in |2C| ]. Note that once a functional form has been chosen for p(a), very different tails can be observed 
depending on the range of variation for a (compare the circles and squares in Fig. ^). This feature might question 
the relevance of the exponent B as an intrinsic quantity for granular gases in steady states. 

An interesting issue concerns the fact that no violation of molecular chaos has been observed: MD and DSMC results 
are in remarkable agreement (even with a packing fraction as high as 40% in MD). This is in contrast with the situation 
of free cooling |n| but also with MD results on heated inelastic hard spheres J3^]. The microscopic precollisional 
velocity correlations driven by the standard inelastic collision rule with a constant restitution coefficient |33|l seem 
significantly reduced within the present random a model. A thorough investigation of short scale velocity correlations 
would require the computation of various precollisional averages involving moments of the relative velocities, and 
has not been performed. Our results however suggest that the dynamical correlations inducing recollisions |j4j and 
responsible for the violation of molecular chaos may not be a generic feature of driven granular gases exhibiting a non 
equilibrium stationary state. 

More refined models could introduce correlations between the effective restitution coefficient and the relative veloc- 
ities of colliding pairs. This feature, neglected here, seems difficult to quantify from first principles, but might affect 
the high energy tail or induce precollisional velocity correlations. It would however be very interesting to be able to 
link a realistic energy injection mechanism with a precise distribution of restitution coefficients. 

Finally, a hydrodynamic study of the present random a model, in which the conservation of the energy is valid on 
average only, is left for future investigations. 



In this appendix, we consider particular distributions p(a) (with a 2 = 1), and give the corresponding formulae for 
a 2 . 



APPENDIX A: FOURTH CUMULANT OF THE VELOCITY DISTRIBUTION 



• Trimodal p(a) = §[<5(a - \/TTj) + 6 (a - vT^y)] + (1 - b)6(a 



1) 



I667 2 



(Al) 



a 2 = 



8&( % /TT7 + x/l - 7) + 32 - 166 - 15bj 2 
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• Bimodal p(a) — \{8{a — + 7) + 5(a — y/l — 7)) (particular case of the trimodal distribution, with b = 1) 

!67 2 

«2 = s f ; (A2) 

16- 1572 + 8(^/1^7 + vT^t) 

Example: for 7 = 0.5, a 2 w 0.1444 

• Flat distribution for a 2 , between 1 — 7 and 1 + 7 

= iey fA , 

^ 16 ((1 + 7) 3 / 2 -(1 - 7) 3 / 2 ) +48 7 - 15 7 3 

Example: for 7 = 0.5, a 2 w 0.0436; for 7 = 1, a 2 » 0.2045. 

• Flat distribution for a, between 7 and 1 + 7 



a 2 = 1 imposes 7 = (y/TT/3 — l)/2, and yields a 2 ~ 0.1868. 

APPENDIX B: GAUSSIAN IFF p = <5(a - 1) 

We show in this appendix that only the trivial elastic distribution p(a) — 8 (a — 1) can lead to Gaussian velocity 
statistics. 

Since we assume a 2 = 1, p,2 = (see section III-A). Therefore, equation (|6|) reduces to 

dap(a)I(f, /) = . (Bl) 

Let us assume that /(c) is Gaussian. The equations for vj* and v 2 * lead to 

a" 2 -l 



/(cD/(c 2 *) = /(ci)/(c 2 ) exp =— (cia ■ &)' ) ■ (B2) 



Equation (Bl) is then recast, by carrying out the integration over a for the term /(ci)/(c 2 ), and simplifying by /(ci), 
into 

Idld J dc 2 c 12 /(c 2 ) = J daf^- J dc2C 1 2f{c 2 )Jd (B3) 

where 

,tt/2 

7^d = 7d / rff sin d ~ 2 6» cos 6> 

J-tt/2 

J d = y dercostfexp ^— — — - — -(c 12 cos #) 2 

/7r/2 / — 2 1 \ 

d0sin d-2 0cos0exp ( Cl2 cos6») 2 J . (B4) 

is the angle between <r and Ci 2 , and jd is a geometrical factor corresponding to the integration over the remaining 
angles. By expanding the exponential, we thus obtain the relation valid for any c±: 



= Id J dc 2 ci 2 f(c 2 ) - J da 



a 2 



/ dc2/(ca) V " cl p 2 +1 / d0 sin^ 2 cos 2p+i (B5 ) 

■/ p=0 V 2 / J-TT/2 

Since this is valid for any Ci, each term of the expansion in powers of c\ 2 must be zero. For p = we obtain 
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a~ 2 = 1 , (B6) 

and for p > 1: 

Q!- 2 (1-Q!- 2 )P = 0. (B7) 

A straightforward recurrence yields a~ 2 P = 1 for any p > 0, and thus 

p(a) = 5(a - 1) (B8) 



p(a) 


a 


a 4 


theoretical ai 


measured ai 


bimodal 1 


0.966 


5/4 


0.144 


0.13 


bimodal 2 


0.999 


1.00666 


3.3 10" 3 


3.3 10~ 3 


flat a 


0.957 


1.31111 


0.187 


0.162 


flat a 2 G [0; 2] 


0.943 


4/3 


0.2 


0.178 


trimodal 1 


0.943 


4/3 


0.2 


0.178 


flat a 2 € [0.5; 1.5] 


0.989 


1.0833 


4.4 10~ 2 


4.2 10~ 2 


trimodal 2 


0.989 


1.0833 


4.4 10~ 2 


4.2 10~ 2 


Gaussian, s = 0.1 


0.986 


1.09997 


5.3 10~ 2 


0.051 


Gaussian, s — 0.2 


0.968 


1.199 


0.11 


0.102 



TABLE I. bimodal l:±(8{a-- s /T/2)+5(a-- s /3/2) ; bimodal 2:±(<S(a-1.04)+<5(a-0.958332) ; flat a: a £ [0.457427; 1.457427] 
trimodal: p(a) = §(<5(a - ^/T+r) + 5(a - V^^)) + (1 - b)5{a - 1); trimodal 1: 7 w 0.835254, b w 0.47779 ; trimodal 2: 
7 « 0.39058, 6 « 0.546248; Gaussian: rfto(a 2 ) oc 9(a 2 ) cxp(-(a 2 - l) 2 /(2s 2 )) 
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